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This technical report presents an algorithm and code to solve for two similar 
classes of problems. Given a list of points, find the equations for the best line 
that passes a minimum distance from each point. Given a list of lines, find 
the point that is a minimum distance from each line. Singular value decom- 
position is used to solve for the inconsistant matrix. Numerical robustness is 
tested for alternate solutions to the line problem, and the best solution form 
is returned. This avoids singularities in best fit lines that are oriented along 
certain directions. This algorithm is implemented in M C" , and has applications 
in locating points in tomographic slice stacks, and finding trajectories of high 
energy particles passing through arrays of detectors. 

The source code is presented in the appendix. 

• Find the line that has the minimum RMS 1 distance from a list of points 
denoted (xi , y, , 2,) , i = 1, 2, 3, . . . 

• Find that point that has the mimimum RMS distance from a given list of 
lines. 

*This work was supported by an equipment loan from Silicon Graphics, Inc. D. Karron is 
with Department of Applied Science at New York University and the Department of Surgery 
at New York University Medical Center, 560 First Avenue, New York, New York, 10016. 

1 Root Mean Squared 
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Abstract 



An algorithm is presented to solve for two similar classes of problems. 
Given a list of points, find the equations for the best line that passes a 
minimum distance from each point. Given a list of lines, find the point 
that is a minimum distance from each line. Singular value decomposition 
is used to solve for the in consist ant matrix. Numerical robustness is tested 
for alternate solutions to the line problem, and the best solution form 
is returned. This avoids singularities in best fit lines that are oriented 
along certain directions. This algorithm is implemented in "C", and has 
applications in locating points in tomographic slice stacks, and finding 
trajectories of high energy particles passing through arrays of detectors. 
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It is crucial to represent each line as the intersection of two planes. The 
solution of the second problem then follows simply from the results of the first 
problem. Each line is given as a pair of intersecting plane equations : 

Azi + B yi + C Zi | ,ss 1 

D x { ; + E yi ; + F Zi = 1 

where t = 1, 2, 3, . . ., n 

The first problem can be formulated as 



[A] [X] = [B] 



where 



[A] = 



*o Vo *o 

*i yi *i 

*2 V2 Z2 



. Xn y n Zn . 



[X} = 



A D 
B E 
C F 



[B] = 



1 
1 
1 

• 

1 1 



Note that the data matrix [A] has, in general, many more rows than columns, 
so that this system is usually over determined. The SVD (Singular Value De- 
composition) solution is thus one of many possible solutions and can be shown 
to be the one with the required RMS error properties. 

The SVD technique is based on the observation that any matrix [A] whose 
number of rows N is greater or equal to its number of columns N , can be 
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written as the product of an /V x M c lumn orthogonal matrix [U], an M x M 
diagonal matrix [W], with positive orj zero elements, and the transpose of an 
M x M orthogonal matrix [V] [1], [2], 3]; that is, 



\ A ]=[dl[W][V]T. 

The "standard" SVD solution, which can be verified by back substitution, 
is given as: 



I*J=[V][WI-1[U]TIB] 

The problem with this solution is that the pairs of planes that result can 
be close to parallel, leading to severe numerical instabilities if these values are 
used to solve the second problem. To minimize these numerical instabilities we 
constrain the solution planes to be perpendicular. Any calculation with these 
lines will suffer from numerical instability if the rows of the line matrix are not 
"maximally" linearly independent. 

Intuatively, think of the line as two intersecting planes. Think of each plane 
as minimally defined as three points. Now at first glance it would appear that 
6 points are required to define two planes. Remember that the planes can have 
points in common, and that we can take two triangles for two planes that share 
two points. So we can minimally solve for two planes in space with 4 points. 
If we now constrain our planes to be perpendicular, then we do not need one 
point. Now if we also allow our perpendicular planes to rotate about the line, 
then we do not need the remaining third point to constrain the perpendicular 
planes in any orientation with respect to our line. Therefore, we can use two 
points to define a line, but the orientations of the planes are undefined, but 
perpendicular to each other. 

From the above argument, it is clear that there are multiple equivalent pairs 
of planar equations for any line. Perpendicularity means that any one coefficient 
from the first plane {A, B, C) and any one coefficient not in the same column 
{D, E, F} can be 0. This is equivalent to projecting the line onto any two of the 
x — y>y — z < or z — x planes. This reduces the 2 x n overdetermined equations in 
6 unknowns to a three systems of n overdetermined equations in two unknowns. 
However, if any of the x; ,y; , or z; values for our lines approach zero, then 
one pair of equations becomes singular. To avoid singularity in any particular 
direction, we must determine which three pairs of equations to use. Since each 
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{ 



b[i]=1.0; /* point normal plane constant */ 



forCj=la<2u++) 
{ 

u[i][j]=a[i][j]; /* don't touch matrix a */ 
} 



Print_Vector(" cons t ant s " , 1 ,b, 1 ,m); 
#cndif DEBUG 

8vdcmp(u,m,n,w,v); /* decompose u into u,v, and w */ 

wmax=0.0; /* edit diagnonal which is vector v */ 
for(i=l;i<m;i++) 

{ 

if(w[i] > wmax) 



} 

} 

#ifdef DEBUG 
printf("wmax=y.f \n" ,wmax); 
#endif DEBUG 

wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 
{ 

w[i]=0.0; 
} 



} 



#ifdcf DEBUG 
I* Lets see the intput */ 
Print-Matrix("- ' 



".l.a.l.m.l.n); 




} 



#ifdef DEBUG 
printf("wm» 
#cndif DEBUG 



n",wmin); 



8vbksb(u,w > v,m,n > b,x); /* svd back substitute to solve for best x +f 
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#ifdef DEBUG 

Print_Matrix( M lhs",l > u,l,m,l l n); 
Print_Vector( M d« ,l,w,l,m); 
Print_Matrix( M rhs M , 1 , v, 1 ,m, 1 ,n) ; 
Print.Vector( M solut» ,l,x,l,n); 
#endif DEBUG 

nrdot(a,l,m,l,n,x,c); 

#ifdef DEBUG 

Print_Vector("test solut» ) , ',l,c,l,m); 
#endif DEBUG 

condition=0.0; f* this is not really condition, just an figure of merit 
* for the absolute error, as absolute difference from 1.0 

*/ 

for(i=l;i<m;i++) 
{ 

if((c[i]-1.0)<0.0) 

condition — = (c[i]— 1.0); 
else 

condition •+•= (c[i]— 1.0); 

} 

#ifdef DEBUG 

printf("cond" n=y.l\n", condition); 
#endif DEBUG 

free_matrix(u, 1 ,m, 1 ,n) ; 
free_matrix(v,l,m,l,n); 
free_vector(b, 1 ,m); 
free_vector(c,l,m); 
free_vector(w,l,n); 

return condition; 

} 

void NueNogginKnocker(int m.float **plist,float lcoeff{2][3]) 
{ 

A 

* Given a list of m point triples (x,y,z) (which actually should be declared 

* (*P)$] out * Ba< * trouble making that work) return the point normal line 
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* coefficients (where the point constant is taken as 1.0, always!) 

* in lcoeff for two point normal planes that describe the line that 

* best fits the list of points in three dimensions. 

* 

register int i; 
float **a,*x[4]; 
double c[4]; 
int cmaxi; 
double cmaxf; 
static int n=2; 

#ifdef DEBUG 
Debug=2; 
#endif DEBUG 

#ifdef DEBUG 

for(i=0;i<m;i++) /* Recite the incomming point list just to be certain. */ 
{ 

printf("IueNogg- nocker Po» , V.i , %f \n" ,i, 

*(plist[i]), *(pli8t[i]+l), *(plistp]+2))j 

} 

#endif DEBUG 

if( m < n ) 
{ 

fprintf(stderr,"too lew po« y.d<y.d to do anyth" M<H\n",m,n); 
return; 

} 

a=matrix(l,m,l,n); 

x[l]=vector(l,n); /* iiang three solution vectors off x */ 

x[2]=vector(l,n); 

x[3]=vector(l,n); 

A first try */ 
for(i=0;i<m;i++) 

{ 

a[i+l][l]= *(plist[i]); /♦ x */ 
a[i+l][2]= *(plist[i]+l);A y */ 
} 

#ifdef DEBUG 
printf("l« try\n"); 
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#endif DEBUG 

c [1] = RunCalculation(a,m,x[l] ) ; 

/* second try */ 

for(i=0;i<m;i++) 
{ 

a[i+l][l]= *(plist[i]); /***/ 
a[i+l][2]= *(plist[i]+2);A z */ 
} 

#ifdef DEBUG 
printf(" second try\n"); 
#endif DEBUG 

c[2]=RunCalculation(a,m,x[2]); 

/* third try */ 

for(i=0;i<m;i++) 
{ 

a[i+l][l]= *(plist[i]+l); /♦ y */ 
a[i+l][2]= *(plist[i]+2); /* z */ 
} 

#ifdef DEBUG 
printf( M th« try\n M ); 
#endif DEBUG 

c[3]=RunCalculation(a,m,x[3]); 

cmaxf=0.0; /* which is the worst of the three ?, keep the two best ! */ 
for(i=l;i<3;i++) 

{ 

if(cmaxf<c[i]) 
{ 

cmaxf==c[i]; 
cmaxi=i; 

> 

} 

#ifdef DEBUG 

printf( M cmax= , /,l , cmax«=y.d\n",cmaxf,cmaxi); 
#endif DEBUG 
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switch(cmaxi) /* wire the best two solutions 

* back into the solution for the line 

V 

{ 

case 3: 

lcoeff[0][0]=x[l][lj; 
lcoeff[0][l]=x[l][2]; 
lcoeff[0][2]=0.0; 

lcoeff[l][0]=x[2][l]; 
lcoeff[l)[l]=0.0; 
lcoeff[l][2]=x[2][2]; 
break; 

case 2: 

lcoeff[0][0]=x[l][l]; 
lcoefF[0][l]=x[l][2]; 
lcoeff[0][2]=0.0; 

lcoeff[l][0]=0.0; 
lcoeff[l][lJ=x[3][l]; 
lcoeff[l][2]=x[3][2]; 
break; 

case 1: 

lcoeff[0][0]=x[2][l]; 

lcoeff[0][l]=0.0; 

lcoeff[0][2]=x[2][2]; 

lcoeff[l][0]=0.0; 
lcoefitl][l]=x[3][l]; 
lcoeff[l][2]=x[3][2]; 
break; 

} 

#ifdef DEBUG 

printf("lcoef f [0] = It ,V.i ,y.l\n",lcoefftO][0],lcoerTtO][l],lcoerTlO][2]); 
printf("lcoef f [1] = */.f ,%1 ,y.l\n",lcoefill][0],lcoerT[l][l],lcoerTll][2]); 
#endif DEBUG 

free_matrix(a, 1 ,m, 1 ,n) ; 
free_vector(x[l] , 1 ,n) ; 
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free_vector(x[2] , 1 ,n); 
frec.vector(x[3] , 1 ,n); 

} 

/****************************************^ 

#ifdef OBSOLETE 

I* this version does not check for best two solutions out of three, 
* and has a bad numerical instability in certain directions. 

v 

void NogginKnocker(int m,float **plist,float lcoeff[2][3]) 
{/* project onto two planes, first pass= x,y */ 
/* second pass, x,z */ 
register int ij,k; 

float **u,*w,**v,**a,*b,*c,*x,**cvm; 
float wmaXjWmin; 
int n=2; 

float al,a2,bl,c2; 

#ifdef DEBUG 
Debug=0; 
#endif DEBUG 

forQ=0J<ma++) 
{ 

printf("Po» id ,id ,Xi\n" j, 

♦(plistDl), *(pli8t[i]+l), *(plist[j]+2)); 

} 

if(m < n ) 
{ 

fprintf(stderr,"Umlout M<I\n"); 
return; 

} 

a=matrix( 1 ,m, 1 ,n) ; 
b=vector(l,m); 
c=vector(l,m); 
u=matrix( 1 ,ra, 1 , n ) ; 
v=matrix( 1 ,m,l ,n) ; 
c vm=matrix( 1 ,m, 1 ,n) ; 
w=vector(l,m); 
x= vector (l,n); 

for(i=0;i<m;i++) 
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{ 

b[i+l]=1.0; 

u[i+l][l]=a[i+l][l]= *(plist[i]); /* x */ 
u[i+l][2]=a[i+l][2]= *(plist[i]+l);A y */ 
} 

Print_Matrix("F- Pass. ■< .l.a.l.m.l.n); 

Print_Vector("F r < Pass, constaats" 1 l,b,l,m); 

svdcmp(u,m,n,w,v); 

wmax=0.0; 

for(i=l;i<m;i++) 

{ 

if(w[i] > wmax) 
{ 

wmax=w[i]; 
} 

} 

#ifdef DEBUG 

printf("F" Pass wmax=y.f\n",wmax); 
#endif DEBUG 
wmin=wmax*(1.0e-6); 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 
{ 

w[i]=0.0; 
} 

} 

8vbksb(u,w > v,m > n,b,x); 
#ifdef DEBUG 

Print_Matrix( " lhs " , 1 ,u , 1 ,m, 1 ,n) ; 
Print_Vector("d- .l.w.l.m); 
Print_Matrix("rlis' , l l > v,l,m,l,n); 
Print.VectorC'F- pass solut" ,l,x,l,n); 
nrdot(a,l,m,l,n,x,c); 

Print_Vector("F« pass test solut- ?) l, ,l,c,l^n); 
#endif DEBUG 

al=x[l]; 
bl=x[2]; 

/* second pass */ 
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for(i=0;i<m;i++) 
{ 

u[i+l][l]=a[i+l][l]= *(pliflt[i]); A x V 
u[i+l][2]=a[i+l][2]= *(plist[i]+2);A z */ 

} 

Print-Matrix("SecondPas8, « ,l,a,l,m,l,n); 

Print_Vector("Second Pass, constants M I l,b 1 l,m); 

svdcmp^.ra.n.w.v); 

wmax=0.0; 

for(i=l;i<m;i++) 

{ 

if(w[i] > wmax) 
{ 

wmax=w[i]; 
} 

} 

#ifdef DEBUG 

printf("Secpmd Pass wmax=y,l\n",wmax); 
#endif DEBUG 
wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 
{ 

if(w[i] < wmin) 

{ 

w[i]=0.0; 
} 

} 

8vbk8b(u l w,v J m,n,b,x); 
#ifdcf DEBUG 

Print_Matrix("llis", l,u f l,m,l,n); 
Print_Vector("d» f l,w,l,m); 
Print_Matrix( ,, rhs ,, ,l l v,l,m > l l n); 
Print_Vector("Second pass solute ,l,x,l,n); 
nrdot(a,l,m l l,n,x,c); 

Print.VcctorC'Second pass test solut- W) m J.# t ljn); 
#endif DEBUG 



a2=x[l]; 
c2=x[2]; 

free_matrix(a, 1 ,m, 1 ,n) ; 
free-vcctor(b,l,m); 
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free_vector(c, 1 ,m); 
free_matrix(u, 1 ,m, 1 ,n) ; 
free_matrix( v, 1 ,m, 1 ,n) ; 
free.vector (w, 1 ,n) ; 
free_vector(x,l,n); 
free_matrix(cvm, 1 ,m, 1 ,n) ; 

#ifdef DEBUG 
Debug=0; 
#endif DEBUG 

lcoeff[0][0]=al; 
lcoeff[0][l]=bl; 
lcoeff[0][2]=0.0; 
lcoeff[l][0]=a2; 
lcoeff[l][l]=0.0; 
lcoeff[l][2]=c2; 

} 

#endif OBSOLETE 
/*****************************^^ 

void NitherNoid(int mm,float **llist,float point [3], int debug_me) 

{ 

/* Find the best point to solve a mm long list of lines attached to Hist. 

* Actually, the solution is defined the same as the best point to solve 

* for a list of lines consisting of pairs of planes. 

register int ij,k; 

float **u,*w,**v,**a,*b,*c,*x,**cvm; 
float wmax.wmin; 

int n=3; /* solution is in 3 dimensions */ 

int m=mm*2; /* each line consists of two planes */ 

#ifdef DEBUG 
Debug=l; 
#endif DEBUG 



#ifdef DEBUG 
if(debug_me) 
for(j=0a<mm-j++) 
{ 

printf( M I« :L« e[ # /.d] : , /.+e,y.+e, , /.+e\n\t , /.+e. , /.+e, , /.+e\n" j, 

*(llist[j]+0),*(llist[j]+l) l *(mst[i]+2), 
*(Hi8t[j]+3),*(lli8t[j]+4),*(lli8t[j]+5)); 

} 
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#endif DEBUG 

if(m < n ) 
{ 

fprintf(stderr, M Too few plane pa" to solve lor M<N\n M ); 
return; 

} 

a=matrix(l,m,l l n); 
b=vector(l,m); 
c=vector(l,m); 
u=matrbc( 1 ,m, 1 ,n) ; 
v=matrix( 1 ,m, 1 ,n) ; 
c vm=matrix( 1 ,m,l ,n) ; 
w=vector(l,m); 
x= vector (l,n); 

for(i=0;i<mm;i++) 

{ 

b[i*2+l]=b[i*2+2]=1.0; /* point normal plane constant */ 

u[i*2+l][l]=a[i*2+l][l]= *(llist[i]+0); 
u[i*2+ 1] [2] =a[i*2+ 1] [2] = * (Hist [i] + 1 ) ; 
u[i*2+l][3]=a[i*2+l][3]= *(llist[i]+2); 

u[i*2+2][l]=a[i*2+2][l]= *(llist[i]+3); 
u[i*2+2][2]=a[i*2+2][2]= *(llist[i]+4); 
u[i*2+2][3j=a[i*2+2][3]= *(llist[i]+5); 

} 

#ifdef DEBUG 

Print_Matrix("H» « put H ,l,a,l,m,l,n); 

Print-Vector (" N« s t ant s " , 1 ,b, 1 ,m) ; 

#endif DEBUG 

8vdcmp(u,m,n,w,v); 

#ifdefSVD_EDIT 
wmax=0.0; 
for(i=l;i<m;i++) 
{ 

if(w[i] > wmax) 

{ 

wmax=w[i]; 
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} 

} 

wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 
{ 

if(w[i] < wmin) 
{ 

w[i]=0.0; 
} 

} 

#endif SVD_EDIT 
8vbksb(u 1 w,v,m ) n l b,x); 
#ifdef DEBUG 

Print_Matrix( H lhs M ,l,u,l,m,l,n); 
Print_Vector("d" ,l,w,l,m); 
Print_Matrix( " rhs M , 1 ,v, 1 , m, 1 ,n) ; 
Print_Vector("H- eolut" ,l,x,l,n); 

nrdot(a,l,m,l l n,x,c); 

Print.Vector("I- test solut- ?)",l,c,l,m); 

#endif DEBUG 

point[0]=x[l]; 
point[l]=x[2]; 
point[2]=x[3]; 

free_matrix(a, 1 ,m, 1 ,n) ; 
free_vector(b,l,m); 
free_vector(c , 1 ,m) ; 
freejnatrix^.l.m.l.n); 
free_matrix(v,l,m,l,n); 
free_vector(w 1 l,n); 
free_vector(x,l,n); 
free_matrix(cvm,l ,m, l,n); 

} 
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